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Abstract. We report part of our recent work on viscous hydrodynamics with consistent phase 
space distribution f(x, p) for freeze out. We develop the gradient expansion formalism based 
on kinetic theory, and with the constraints from the comparison between hydrodynamics and 
kinetic theory, viscous corrections to f(x,p) can be consistently determined order by order. 
Then with the obtained f(x, p), second order viscous hydrodynamical calculations are carried 
out for elliptic flow «2. 



1. Introduction 

Relativistic hydrodynamics is an important theoretical tool in heavy ion collisions. Especially it 
successfully reproduces the observed elliptic flow V2- In the past several years, with respect 
to the updated understanding of collective phenomena in heavy-ion collisions, two major 
progresses have been made in hydrodynamical models. One is the realization of fluctuations 
in initial state [1], which leads to the practically used event-by-event hydrodynamics [2, 3j. And 
correspondingly many new types of anisotropic flow, such as dipole flow t>i[3] and triangular 
flow f3[U [5], are studied. The other one is the inclusion of viscous corrections [6]. 

In practice, although ideal hydrodynamics has been found effective in characterizing the 
collective medium expansion, viscous corrections are not negligible [7] . One particular example 
of viscous corrections to hydrodynamics is the viscous damping of elliptic flow [HI \9\ llOj. 

Many attempts have been made in the investigation of the viscous corrections to 
hydrodynamics [HI EJ [TO]. Since first order viscous hydrodynamics suffers from the causality 
problem(see for example [11]), most of the established viscous hydrodynamical models have 
contained second order viscous terms in the equations of motion, such as Israel-Stewart theory. 
Recently, in an important paper, Baier and his collaborators(BRSSS) developed second order 
viscous hydrodynamical theory with conformal symmetry assumption[12j. Hydrodynamical 
simulation to the medium expansion in heavy ion collisions is achieved by solving hydrodynamical 
equations of motion, with respect to the specified initial conditions. Observables are obtained 
then by freeze out process. Due to the correspondence between hydrodynamics and kinetic 
theory, freeze out process also needs viscous corrections, and the corrections have to be 
constrained so that they are consistent with hydrodynamics itself. In particular, in Cooper- 
Fryer formula, this is reflected in the determination of the phase-space distribution function 
f(x,p). However, to the knowledge of authors, this constraints have not been systematically 
discussed in the literature through second order, and this will be the main subject of this paper. 
For more details, refer to |13j . 

*In collaboration with Derek Teaney. 



The paper is organized as follows. The theoretical formalism is constructed by reviewing 
viscous effects in kinetic theory and hydrodynamics in section [2j In section [3] the determination 
of phase-space distribution function is discussed. Hydrodynamical simulations with consistent 
form of f(x, p) is carried out, and elliptic flow v% is calculated at RHIC energy, in section |4j 
Our conclusions are summarized in section [H 

2. Kinetic Theory and Hydrodynamics 

A liquid system out of equilibrium can be approached through either transport theory or 
hydrodynamics, provided that the non-equilibrium part can be seen as perturbations. As 
a result, in order to sustain the continuity around freeze out, there must be a determined 
correspondence between kinetic theory and hydrodynamics. This correspondence originates 
the constraints on distribution function f(x,p). In terms of the dependence on transverse 
momentum pp, one can generally fix 5f by the moment method [H], in which 5f is decomposed 
into moment expansion and cut at a certain order. The widely accepted form, for example, 

_ n p (l ± n p ) 

2(e + V)T 2 p p { ' 

is obtained by a fourteen-momentum method. However there are two flaws in this form for 
realistic simulations. First, in actual simulations tt^ is taken from second order hydrodynamical 
calculations, while 5f n is computed using a first order approximation. Second, and more 
importantly, the pt dependence in 6f has been investigated in [15], and found to vary based on 
microscopic dynamics. Although the quadratic ansatz (i.e. Eq. ([Tj)) is appropriate for collisional 
energy loss, it fails in the case where radiative energy loss dominates. 

In kinetic theory, the viscous corrections Sf(x, p) corresponds to the viscous corrections to 
hydrodynamics, i.e., stress tensor tt^ v , so it also depends on gradient expansion. In this way, 
following Chapman-Enskog method pH] where expansion in gradients is used to find asymptotic 
form of the solution to Boltzmann equation, 5f(x, p) can be determined order by order. When 
relaxation time approximation can be taken into account, form of 5f is analytically solvable for 
the specified viscous hydrodynamics. In this paper, we take relaxation time approximation and 
determine 5f(x,p) through second order in gradients, with respect to BRSSS hydrodynamics. 

The matrix convention we use throughout the paper is (—,+,+,-1-). p^ is used for four- 
momentum. And for the sake of convenience, we introduce a tilde on a quantity from time to 
time to indicate it as dimensionless. In the derivation of hydrodynamics, tensor index is always 
split into temporal and spatial parts. Since flow four-velocity u^ 1 = j(l,v) is purely temporal in 
the local rest frame(LRF), the decomposition is then realized based on and the projection 
operator = u u u u + g^ u '. For derivative operator we have = u^D + V 4 , with D = u^d^ 
and = A^du representing pure time and space derivatives in LRF respectively. A further 
decomposition according to rotation group is also considered comparing to the existed tensor 
form in hydrodynamics, so the indices of a tensor in brackets (. . .) stand for being symmetric, 
(projection operator) traceless and transverse, 

jyi...(...a...p...)...v _ j^i...{...p...a...)...v 

A»-(- a -)~»A a p = 
More details of the decomposition can be found in [13] . 

2.1. Hydrodynamics 

In Chapman-Enskog expansion, the so-called solubility condition[14J relates spatial gradients 
and temporal derivatives. This is actually equivalent to hydrodynamical equations of motion 



(2) 
(3) 
(4) 



df Ji T fJ,L ' = 0. In the current section we will derive the viscous corrections through second order 
to BRSSS hydrodynamical equations of motion. 

The energy-momentum tensor of viscous hydrodynamics reads 

= eu»u u + VA^ U + vr^ , (5) 

where ■n^ v is the viscous correction to stress tensor. In this paper, the bulk viscosity Q = 
and baryon chemical potential [ib = 0. Up to second order viscous corrections, BRSSS [16 
determined that the possible forms of the gradient expansion in the stress tensor in a conformal 
liquid are 



+Ai ( a\a vX ) +A 2 ( a\W x ) +A 3 ( n\W x 
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where a^ u = 2(W^u v ) and the vorticity tensor is defined as 

W = I A ^A^ (V a up - Vpu Q ) . (7) 

(r n , Ai, A2, A3) are corresponding second order transport coefficients, d = 4 is the number of 
space-time dimensions jf] Then equations of motion according to Eq. ^ can be formulated as 

De = - (e + V)V ■ u + -v^c^ + 0(V 3 ) , (8) 

Du, = - + U ^ A - - A - V ^ + Q(V3) . (9) 

M e+P e+P e+P V 7 V ; 

Note that the temporal derivatives of the flow velocity and the energy density e are written 
in terms of spatial gradients through second order, while higher order terms are neglected. 

2.2. Kinetics 

To determine the viscous corrections to the distribution function, the strategy is to solve the 
kinetic equations in a relaxation time approximation order by order in the gradient expansion: 

f(x, p) = n p (x, p) + 5 f(x, p) = n p (x) + df {1) + 5f {2) + 0(V 3 ) , (10) 

where n p (x) is equilibrium distribution, which depends on the space-time coordinates through 
the dimensionless combination p ■ u/T = p ■ u. In classical limit, n p (x) = exp[p • u(x)/T(x)}. In 
a relaxation time approximation the Boltzmann equation reads 

P^f p (x) = - ^ [/ - n{p ■ u*(x)/T m (x))] , (11) 

where the dimensionless coefficient C p is related to the canonical momentum dependent 
relaxation time tr, 

C p = —T 2 TR{ ~ P - U) . (12) 
p ■ u 

fin this paper, we will ignore the possible effects of higher space-time dimensions, so we write explicitly 
wherever d = 4. Also we ignore all curved space-time effects discussed in [To] . 



Following the convention used in |15j . the canonical relaxation time is proportional to (— p-u) 1 ^ 01 , 
where a lies generally between the quadratic ansatz limit(a = 0) and linear ansatz limit (a = 1). 
We have 

C p = c (-p-u)- a , (13) 

with Co to be fixed by shear viscosity (see below). The parameters u*^{x) and T*(x) which 
appear in the Eq. (11) are equal at leading order to the temperature and flow velocities (T and 

(14) 



u^). They satisfy the Landau matching conditions, i.e. , 
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So there exist the expansion 



p^5f(u*,T*) < = 0, 
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T*{X) =T(X)+5T* , 
u*"(X) =u^{X) + 6u*» . 
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with 5T* = Xm =1 5T( n j and 5u*^ = Y2 n =i ^ u (n) i n dicating corrections from all higher order 
gradients. Expanding n* (with an obvious notation) we find correspondingly, 



n p + 5n* and <5n* 



p ■ 5u* 
T 



5% 

P rp2 



(17) 



where here and in the following a prime stands for the derivative on p ■ u. Jn* with higher 
order gradients is important for the determination of viscous correction 5f, as will be shown 
in the next section. Eq. (14), together with the correspondence between kinetic theory and 
hydrodynamics , 

' /,/; -p^5f(u*,T*) = 7r^, (18) 



(2^) 3 £p 

constrains the solution of Boltzmann equation /(x,p). 
3. Determination of Sf 



We substitute the expansion Eq. (10) into the kinetic equation Eq. (11) and equate orders. In 
doing so we use the hydrodynamic equations of motion to write time derivatives of T{x) and 
u(x) in terms of spatial gradients of these fields, see Eq. ([8]). Indeed, with the help of these 
equations of motion, up to first order in gradient, 



p ^d tl n p = +n' p ^ r a^ + 0(V) 



(19) 



The terms linearly proportional to the o^ v is ultimately responsible for shear viscosity. 

Collecting all the terms of first order in gradient for 5f we obtain the preliminary first order 
solution 



C P n p 
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p(i)- 



(20) 



But Eq. (20) is not completed until the constant coefficient Co, the form of Su^ and ST^ 
are determined through constraint conditions Eq. (14) and Eq. (18). Obviously there exist the 
trivial solution 5u^ = and ST(i) = to Eq. (14), and then from Eq. (1 1 Sh we find, 



(21) 



B\ is a dimensionless constant from the integral, 



5i(a) = v 



' /3 V(p-«)HA.j 



(2tt)3^ 



(22) 



It only depends on the underlying microscopic dynamics in terms of a. We will encounter three 
more similar constants, 
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Co is then the solution to Eq. (21), and it is interesting to know that for a massless gas and in 
the quadratic ansatz limit, cq = rj/s. 

We extend the derivation to second order #/(2), then 
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Clearly there are three sources of contributions, those generated from <5/(i) in Boltzmann 
equation, from second order corrections to flow velocity and temperature and from second order 
hydrodynamics equations of motion. The dominant term is from p^d^Sfru, since in Boltzmann 
equation the derivative gives rise to one higher order dependence on pj< simultaneously. To write 
8f(2) explicitly we decompose the resulting tensors into irreducible tensors of the rotation group 
in the LRF ; |T3] provides a few more details. Straightforward then from (26), with somewhat 
tedious algebra, we obtain, 
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The first line in (27) does not contribute to 7r^ u , since all tensor structure from momentum 
integral are orthogonal to these terms in brackets. Besides the already fixed parameter cq, we 
still need to determine second order corrections to flow velocity and temperature, by solving 



Landau-Lifshitz matching condition. And we find that these following solutions satisfy all the 
constraint conditions, up to second order in gradient, 
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As a by-product of the constraint condition Eq. (18), we also recognize that 

~2B 2 



T^B, 



15 
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+ 1 



A 2 



-2r]T w , A 3 



0. 



(30) 



These relations are in the consistent range with those discussed in |17j . In particular, all these 
non-trivial second order transport coefficients are non-free parameters in our formalism, as long 
as r]/s is fixed. This reduces the number of free parameters in hydrodynamical calculations, and 
thus quite significant for the recent study on the extraction of r//s[18j. 

4. Viscous Hydrodynamic Simulation and Elliptic Flow 

The simulation of hydrodynamics in code needs solving hydrodynamical equations of motion. 
For viscous hydrodynamics, for the convenience of computation, the algorithm is introduced by 
approximation. In BRSSS formalism, up to second order in gradients, Eq. ([6]) can be rewritten 

as 



7T 
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such that a linearized equation for tt^ u , which is approximately correct if all higher order 
gradients are neglected, is obtained. Then together with general hydrodynamical equations 
of motion, d^T^ u =0, and appropriately selected equation of state, we have a complete set of 
equations for the unknown variables, (tt^ u , e,u^). There are two options for the equation of state 
in our calculations. For theoretical interest we can consider conformal equation of state(CEOS). 
Although CEOS does not have the accurate description for real medium expansion in heavy ion 
collisions, it has the conformal symmetry as assumed in BRSSS hydrodynamics. The previously 
used lattice equation of state(LEOS) by Romatschke and Luzum[9J is also tested. The details of 
the algorithm can be found in |13j . For the initial state, smooth initial condition is considered. 
As indicated by several groups [31 [T9]. event fluctuations in the calculations of v 2 can be safely 
ignored. A constant temperature freeze out scheme is taken, with Tj Q = 150 MeV. 

The effects of including consistent 5f in freeze out procedure can be best seen in the results 
of elliptic flow. In Fig. [T] the results of a series of calculations are shown. Compare to ideal 
hydrodynamics, viscous hydrodynamics indeed damps the generated v 2 , and a major fraction 
of the damping is from 5f. This is reflected in the difference between the calculated v 2 with 
equilibrium distribution function n p and that with n p + Sf's. The complete and consistent 
calculation with respect to BRSSS hydrodynamics needs also Sft2) correction. However, in 
the quadratic ansatz limit and conformal equation of state the corresponding v 2 (left panel in 
Fig. [TJ quickly becomes unreliable when px goes up to 1.5 GeV. This is due to the dramatic 
increase of viscous corrections at large pt region, and has already been discussed in the first 
order case with <5/,r[l0]. But now that 8f( 2 ) has higher pt{Vt) dependence, this effect is much 
stronger. In another way, we can estimate the magnitude of viscous corrections by calculating 
the spectrum. Then gradient expansion formalism fails when 5f contributes to spectrum no 
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Figure 1. Differential v<i results from various viscous hydrodynamical calculations, b = 6.8 
fin. 
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Figure 2. Spectrum from various viscous hydrodynamical calculations, b = 6.8 fm. Note that 
a factor of 10 has been multiplied in the case of LEOS to make it on the same scale. 



longer perturbatively. As seen in the left panel in Fig. [2j we find similarly for CEOS and 
quadratic ansatz limit, px > 1.5GeV region is beyond the feasibility of our formalism. 

Changing equation of state or the a dependence, however, can help solving the problem to 
some extent. As seen in the two other panels in Fig. [T] and Fig. [2j viscous corrections are smaller. 
On one hand these changes extend our reliable calculation to higher pr region. On the other 
hand, realistic medium in heavy ion collisions has equation of state more close to LEOS, and 
also quadratic ansatz limit does not necessarily reflect the real microscopic dynamics. So as 
expected, second order corrections in 5f lead to extra corrections to calculated observables. 



5. Discussions and Future Work 

We have presented a formalism for the derivation of viscous corrections to phase space 
distribution function at freeze out. With respect to BRSSS hydrodynamics, from this formalism 
we have obtained the consistent form of $f(2) ; an d checked its impact on hydrodynamical 
simulations. &f(2) affects the results as corrections, and the corrections increase with pp. But for 
LEOS or a = 1, the corresponding corrections are as small as perturbations. This is expected 
from the gradient expansion of the formalism we are following. At last, this dependence on a 
verifies the claim in [15] that quadratic ansatz may be questionable in real calculations. 

One important property of higher order viscous corrections in /(x,p), is its higher px 
dependence. This plays a crucial role in the study of anisotropic flow for two aspects. 1. 
This makes the flow results more sensitive to r//s. 2. This makes higher order anisotropic flow 
v n more sensitive. Both of these aspects have positive effects on the extraction of r//s from 
heavy ion collisions. And this will be the main subject of our future work. 
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